Matrix Product States for dynamical simulation of infinite chains 
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We propose a new method for computing the ground state properties and the time evolution of 
infinite chains based on a transverse contraction of the tensor network. The method does not require 
finite size extrapolation and avoids explicit truncation of the bond dimension along the evolution. 
By folding the network in the time direction prior to contraction, time dependent expectation 
values and dynamic correlation functions can be computed after much longer evolution time than 
with any previous method. Moreover, the algorithm we propose can be used for the study of some 
non-invariant infinite chains, including impurity models. 
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Numerical simulation has become a fundamental tool 
to study quantum many-body systems in condensed mat- 
ter physics. Unfortunately, the exponential scaling of the 
dimension of the Hilbert space with system size means 
that brute-force methods are only practical for very small 
system sizes. However, other techniques based on Matrix 
Product States (MPS) [H, i, i, i, Q , such as DMRG 0, 
achieve excellent results in one dimension. 

The success of DMRG methods is based on the fact 
that many interesting physical states, including the 
ground states of many local Hamiltonians, can be well 
approximated by a MPS. For a chain of N d-dimensional 
systems, this has the form 
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Each A\. is a ZJ-dimensional matrix. An important fea- 
ture of MPS is that expectation values of local operators 
can be efficiently computed, allowing them to be used 
variationally. DMRG excels in the computation of static 
properties of finite chains with local Hamiltonians, where 
the required bond dimension grows slowly with the size of 
the syste m M , Q , and of translationally invariant infinite 
chains [olfloIlT However, these methods can break down in 
time-dependent problems far from equilibrium and also 
encounter difficulties in dealing with infinite chains con- 
taining impurities. 



MPS algorithms for non-equilibrium dynamics [ill. |12 



13j, [ij, ll5| work well when the system is close to its 



ground state, but when the system is far from equilib- 
rium the entanglement entropy may grow linearly in time 
and the dimension D required to describe the system will 
grow exponentially [3, [o, , causing these methods to 
break down. Although improved algorithms have been 
developed based on finite propagation speed of correla- 
tions [l9j , all known methods are limited to special cases 
or short times. 

Here we propose an alternative method to compute 



FIG. 1: Standard time evolution with MPS. We start with 
a MPS state, represented here by a chain of circles (tensors) 
connected by summed indices, with open lines for the physical 
spin indices. On this state, a sequence of MPO is applied for 
each step of evolution (1), and the result is truncated (2) to 
the maximal bond dimension D. After iterating the evolution, 
expectation values are computed (3) in the final state. 



both ground states and dynamical quantities for infi- 
nite chains within the MPS formalism. Based on the 
transverse contraction of the tensor network, it allows 
the study of problems not accessible by other methods, 
such as an impurity in an infinite system. It enables 
the calculation of time-dependent expectation values of 
local observables and of few body correlation functions 
at different times, in a new much simpler way. More- 
over, by a folding of the network described below, the 
new method enables dynamical studies that range much 
further in time than any other existing method. 

The standard way of computing dynamical quantities 
with the MPS formalism starts with a state that is (ex- 
actly) described by a MPS H]). Then, some evolution 
operator is applied to it for a given time, making use of a 
Suzuki- Trotter expansion [20| of the total evolution oper- 
ator. Within each discrete time step, the evolution opera- 
tor is broken down into a product of operators. In partic- 
ular, for a nearest neighbor Hamiltonian H — hi^i+i, 
we may write e''"^ « ^-iH^5/2^-iH^s^^iH,s/2 ^ w\iem 

(Ho) contains the hi i^i terms with even (odd) i, so that 
each exponential factor is a product of mutually commut- 
ing local terms. Alternatively, the evolution operator can 
be decomposed as a product of translationally invariant 
Matrix Product Operators (MPO) [2l|. The action of 
one step of evolution on the MPS can be computed by 
applying the corresponding sequence of operators (Fig. [1]) 
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(a)Transverse contraction along space 
direction renders a finite 2D network. 



(b)Transverse 
contraction of 
folded network. 



FIG. 2: Expectation value {0{t)) with the basic transverse 
method (a) and with folding (b), where operators for the same 
time step are grouped together in a double effective operator. 



to yield a MPS with larger bond dimension. This must 
be truncated to keep the best MPS description of the 
evolved state with fixed dimension, D. After repeating 
this procedure for the required number of steps, expec- 
tation values can be calculated in the evolved state. The 
accuracy of the description will however drop exponen- 
tially with the successive truncations. 

Our new method avoids this explicit truncation on the 
bond dimension of the evolved MPS. The basic idea is to 
look at the quantity that we want to compute, say the 
time dependent expectation value of some local operator, 
(\I'(t)|0|^'(t)), as the contraction of a two dimensional 
tensor network, and perform it, not along time, but in 
the direction of space (see Fig. 2(a) I. To construct the 



network, we start from the initial MPS and, for every 
evolution step, apply the proper MPOs. Repeating this 
for the required number of evolution steps, we construct 
the exact evolved MPS (within the Trotter approxima- 
tion), as no truncation is carried out. Finally, we apply 
the local operator O and contract with the Hcrmitian 
conjugate of the evolved state as constructed before. 

The procedure above produces a two dimensional net- 
work, infinite in the spatial direction as the original MPS, 
but finite along the time direction. The expectation value 
we want to compute can now be written as Q , 
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where E{t) — J^i^^i^) ® ^'(^) the transfer matrix 
of the evolved state, Eoit) = J^,j[A'{t) ® {t)]{i\0\j) 
contains the only application of the single-body operator, 
and the bracketed supcrindices on each transfer matrix 
indicate the site of the chain. For a translationally in- 
variant MPO representation of the evolution operator, 
the network retains the invariance [2^ and the transfer 
matrix is the same on every site, except for the single one 
on which O acts [2^. If the largest eigenvalue of E{t), A, 



is non-dcgencrate, E (t) 
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may then substitute the left and right semi-infinite lat- 
tices at both sides of the operator by the left and right 
eigenvectors of E{t) corresponding to the largest eigen- 



value, {L\ and \R), 
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We now specify the algorithm for computing time- 
dependent expectation values in translationally invariant 
infinite chains. The first step is to find the best MPS ap- 
proximation, with given bond dimension D, to the domi- 
nant eigenvectors of E(t). To this end, we repeatedly ap- 
ply the transfer matrix E{t) (already written as a MPO 
along the time direction, see Fig. 2(a) I to the left and to 
the right of an arbitrary initial MPS vector and truncate 
the result to the chosen D, using the technigue for two 
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dimensional tensor networks introduced in 
til convergence is achieved. The procedure yields a MPS 
approximation to the eigenvectors, with the truncation 
taking always place in the space of transverse vectors. 
The second step, computing the numerator and the de- 
nominator in ([2]), can be done very efficiently, as each 
term is a contraction of a MPO acting between a pair of 
MPS. The adaptation of the method to the case of imag- 
inary time evolution is straightforward, so that it is also 
useful for finding ground state properties. In this case, 
the eigenvector calculation is similar to that in transfer 
matrix DMRG algorithms for thermal states [2^ . 

With this approach, we study an infinite chain with an 
impurity. We consider an Ising chain. 
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with gi = 1 Vi 7^ 0, and the impurity represented by a 
different value of the field at site i = 0, go- The system 
is started in a product MPS and imaginary time evo- 
lution is applied for a long time, so that we approach 
the ground state. Then we compute the site dependent 
magnetization, {ai^ ) . Such a calculation cannot be easily 
done with a purely invariant method as iTEBD [lo| , be- 
cause the presence of a singular site will affect a cone of 
tensors as time increases. However, with the transverse 
method, the computation of {L\ and \R) is not modified 
by the presence of the impurity Thus the cost of 

computing the expectation value of a local operator act- 
ing on the position i = in the ground state of this chain 
will be the same as in the translationally invariant case, 
while applying the operator at « 7^ will reduce to the 
contraction of a 2D tensor network of width i + i (Fig. [3]). 

The capabilities of the transverse method regarding 
real time evolution can be further illustrated by the com- 
putation of two-body correlators at different times. If we 
consider two different times, t2 > ti, we may write 

( vi/(0) |ot"l(t2)o["+^l(ii)|*(0)) 

= (vE'(o)|;7(t2,o)tot"ic/(t2,ti)of+''ic/(ti,o)|*(o)) 
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(a)Tensor network 
corresponding to the 
magnetization at distance 
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FIG. 3: Ising chain with magnetic impurity at the origin. 



where E is the transfer matrix resulting from evolution 
until time t2, and Eo-(t) £^re the corresponding MPO con- 
taining the action of each single-body operator [2^ . 

In particular, if both operators act on the same site 
(A = 0), the computation has the same cost as one sin- 
gle expectation value. If A ^ 0, computing ([4]) requires 
instead the contraction of a two dimensional network of 
width A-l-3. This is done by applying one MPO at a time 
and truncating to the closest MPS with the given bond 
D. Since the network is now finite in both directions, 
this last phase of the contraction can be done cither in 
the spatial or in the time direction. 

The success of the transverse approach will depend on 
whether the transfer matrix of the evolved MPS has a 
non-degenerate dominant eigenvector which can be ap- 
proximated by a MPS of reduced dimension. Our imple- 
mentation shows that the procedure achieves comparable 
results to the standard contraction [lo| in a translation- 
ally invariant chain. The transverse method offers the 
advantage of being applicable to dynamical situations in 
which translational symmetry is broken by a small num- 
ber of sites, such as a chain with impurities, or a semi- 
infinite system, but it is also limited to short times. 

However, there is a more efficient representation of the 
entanglement in the transverse eigenvectors. In the MPO 
representing the transfer matrix of the evolved MPS, ten- 
sors that lie at the same distance from the center (occu- 
pied by the physical operator O as in Fig. 2(a) I corre- 
spond to the same time step, coming from a certain term 
and its adjoint in the Trotter decomposition. We can 
group such pairs together in a new MPO by "folding" the 
original MPO (see Fig. |2(b)"| ). The folding operation can 
be understood as performing the equivalent asymmetric 
contraction (^'(t)|0|*(i)) = ($| (^0|*(i))(g)|4'(i))^ where 
4'(t)) is the complex conjugate of the evolved vector and 

\^) = ®k^'i,^i \ikik) is the product of (unnormalized) 
maximally entangled pairs between each site of the chain 
and its conjugate. In our scheme, the ket is now the 
tensor product of two tensor networks corresponding to 
4'(t)} and its conjugate. We may then group together 



each tensor in |^') with the corresponding one in \^), 
and define an effective tensor network of higher bond di- 
mension and physical dimension cP , which can now be 
contracted using again the transverse technique. 

This folded transverse method allows us to explore the 
dynamics until much longer times than any other pro- 
cedure. We may get some physical intuition for this im- 
provement by looking at a single localized excitation that 
propagates freely with velocity v. After time sites xivt 
in the evolved state become entangled. If we look instead 
at the transverse MPS obtained contracting the network 
from the right until x + vt, it is easy to see that all time 
sites are in a product, except for those corresponding to 
the instant t. These sites occupy symmetric positions 
around the center of the network, so that folding groups 
them together in a single site which will be in a product 
state with all the rest. 

As a first benchmark for the new method, we simulate 
the dynamics of states far from equilibrium under the 
Ising Hamiltonian ^ with uniform magnetic field g. The 
initial state |^'o) = ®i-^(|0).i + is evolved with a 
constant Hamiltonian and the results of the transverse 
method with and without folding are compared to the 
exact results (Fig. [5]). For very short times the Trotter 
error dominates in both methods. However, while for 
the transverse procedure (as for iTEBD) truncation error 
becomes soon dominant, and the results deviate abruptly 
from the exact solution, the accuracy of the folded version 
is maintained for much longer times. 

To test the method on a more general problem, we 
repeat the test for a non-integrable Hamiltonian, H = 
— (J2^ <^Wl'^^ + 9<^x + . For this case there are no 
exact results, but we may compare the folded compu- 
tation to the iTEBD simulations with a similar Trotter 
error (Fig. |4]). Again we check that the accuracy of the 
folded procedure for comparable bond dimension reaches 
much longer times. Moreover, remarkably enough, even 
when the results from the folded method start deviat- 
ing (from those to which iTEBD converges for large D), 
they do so in a smooth way, so that, in contrast to other 
procedures, they continue to qualitatively describe the 
evolution for long times. 

This can be seen in a more precise way by looking at 
the truncation error. At a certain time, this can be esti- 
mated by looking at the error in the right eigenvector for 
a given bond dimension with respect to the best eigen- 
vector obtained, i.e. that for the highest D. If we plot 
(Fig. |4]) the bond dimension required to achieve a fixed 
truncation error, we observe that, although the D re- 
quired for a high precision grows exponentially with time, 
with a relatively low bond D < 100 a qualitative descrip- 
tion of the dynamics is reproduced, that lies within 1% 
of the exact solution well beyond times t > 10. 

From the discussion above the transverse method, 
combined with the folding technique, represents a very 
promising tool for the dynamical studies of one dimen- 
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FIG. 4: Magnetization as a function of time. For the Ising 
model (up) results for the transverse method (triangles) are 
compared to the folded version (stars) for D = 60, 120. The 
relative error with respect to the exact result (solid line) is 
shown in the inset. For the non-integrable model (down), 
results with the folded approach for D=60 (red), 120 (blue), 
240 (green) are compared to those of iTEBD (solid lines) for 
increasing values of D. In the inset, the required value of D 
as a function of time, for different levels of accuracy. 



sional systems. The first results show the applicability 
of the method even to non-integrablc systems, allowing 
the simulation of longer evolution times than any other 
technique, and a qualitative description of the dynamics 
until even later. This opens the door for the study of 
physical problems not accessible until now for numerical 
methods, including the dynamics of phase transitions, 
out-of-equilibrium states and thermalization problems. 
The present formalism might also prove very valuable in 
the context of extracting spectral information for quan- 
tum impurity problems, the central problem in dynamical 
mean field theory. The big advantage of our method is 
that we can deal with real frequencies, and no analytic 
continuation from imaginary frequencies is needed as in 
the case of Monte Carlo simulations. In contrast, the 
main limitation would be its exclusive applicability to 
one dimensional systems. Finally, although the method 
has been described for infinite chains, it is easy to adapt 
the technique for the dynamical study of finite systems. 

Note: An independent derivation in |i2()j, in the con- 
text of concatenated tensor network states, lead to a sim- 
ilar network to describe time-evolved states. 
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Using the standard decomposition of the evolution op- 
erator into commuting products, the resulting network 
has translational symmetry with period two, so that the 
argument can be easily adapted substituting E for the 
product of two contiguous transfer matrices, EeEo- 
In the infinite limit we do not need to consider the vector 



terms at the edges, e\ 
different for the finite case. 



and E^\ which would also be 



